####################################################################################################
## State Prohibition Voting Analysis ###############################################################
## Michael Olson and James Snyder ##################################################################
####################################################################################################

###################################################################################################
# Load Assembled Data #############################################################################
###################################################################################################

  pooled <- read.csv("./state_votes.csv")
    pooled$X <- NULL
  
  # drop the counties that we've found change a significant amount in area

  pooled <- pooled[pooled$toobigchange==0,]
    pooled$toobigchange <- NULL
  
  # drop observations with no measure of district preferences
  
  pooled <- pooled[!is.na(pooled$score_mean),]
  
  # drop observations with no vote for prohibition
  
  pooled <- pooled[!is.na(pooled$vote_prohibition),]
  # this choice basically means that the district outcome will be the mean vote choice 
  # of those who cast a vote, rather than the whole district's outcome being marked as
  # missing just because of one person not being there.
  
    # what states are there? we include this in the text
  
      states <- sort(as.character(unique(pooled$long_state)))
      states_string <- paste(paste(states[1:(length(states)-1)],sep=", ",collapse=", "), ", and ", states[length(states)],sep="")
      
      fileConn<-file("./results/included_states.tex")
      writeLines(as.character(states_string), fileConn)
      close(fileConn)
      
  # mark southern states
    
    south <- c("AL","AR","FL","GA","LA","MS","NC","SC","TN","TX","VA")
    
    pooled$south <- 0
    pooled[is.element(pooled$state,south),"south"] <- 1
  
    pooled$Region <- as.factor(ifelse(pooled$south==1,"South","North"))
    
  # aggregate the data for analysis

    mean(pooled$vote_prohibition) # what is the overall support for prohibition?
    
    pooled_ag <- as.data.frame(pooled %>% group_by(state,county,south,score_mean,score_mean_nlo,Region) %>% summarise(vote_prohibition=mean(vote_prohibition),
                                                                                                                      n_votes = dplyr::n()))
    
  # mark the terciles of the district preference score
    
    pooled_ag$terciles <- factor(ifelse(pooled_ag$score_mean<=quantile(pooled_ag$score_mean,0.3333333)[1],1,
                                  ifelse(pooled_ag$score_mean>=quantile(pooled_ag$score_mean,0.6666666)[1],3,2)))
    
  # now match and create a matched dataset
  # we match on the aggregated dataset and then pull out the matched districts for both the aggregated and non-aggregated data
    
    set.seed(53074)
    
    pooled_ag_match <- pooled_ag
    
    south_match_out <- Match(Tr=pooled_ag_match$south,X=pooled_ag_match$score_mean,
                             caliper=0.1,Weight=2,ties=FALSE,
                             distance.tolerance = 1e-20)
  
    included_dists <- pooled_ag_match$county[c(south_match_out$index.control,south_match_out$index.treated)]
    
    list_binder <- list()  # for the aggregated data
    list_binder2 <- list()  # for the non-aggregated data
  
    for(j in 1:length(included_dists)){
      
      list_binder[[j]] <- pooled_ag[pooled_ag$county==included_dists[j],]
      list_binder2[[j]] <- pooled[pooled$county==included_dists[j],]
      
    }
  
    matched_pooled <- do.call(rbind,list_binder)
    matched_pooled$temp_outcome <- matched_pooled$vote_prohibition
    
    matched_pooled_noag <- do.call(rbind,list_binder2)
    matched_pooled_noag$temp_outcome <- matched_pooled_noag$vote_prohibition
    
    pooled_ag$temp_outcome <- pooled_ag$vote_prohibition
  
# match separately for the "no local option" measure
  
    pooled_ag_match_nlo <- pooled_ag[!is.na(pooled_ag$score_mean_nlo),]
    
    south_match_out_nlo <- Match(Tr=pooled_ag_match_nlo$south,X=pooled_ag_match_nlo$score_mean_nlo,
                             caliper=0.1,Weight=2,ties=FALSE,
                             distance.tolerance = 1e-20)
    
    included_dists_nlo <- pooled_ag_match_nlo$county[c(south_match_out_nlo$index.control,south_match_out_nlo$index.treated)]
    
    list_binder_nlo <- list()
    
    for(j in 1:length(included_dists_nlo)){
      
      list_binder_nlo[[j]] <- pooled_ag_match_nlo[pooled_ag_match_nlo$county==included_dists_nlo[j],]
  
    }
    
    matched_pooled_nlo <- do.call(rbind,list_binder_nlo)
    matched_pooled_nlo$temp_outcome <- matched_pooled_nlo$vote_prohibition
    
# make row ids for "clustering" (how we're going to do regular robust SEs) for all the different datasets

    # this is basically a "hack"--clustering by rowids produces "clustered" standard errors that 
    # are identical to heteroskedasticity robust standard errors, but the felm() function automatically
    # attaches clustered ses to the object (so that they get reported in tables without any extra lines of code)
    # which isn't automatic for robust standard errors.
    
    pooled$rowid <- 1:nrow(pooled)
    pooled_ag$rowid <- 1:nrow(pooled_ag)
    matched_pooled$rowid <- 1:nrow(matched_pooled)
    matched_pooled_noag$rowid <- 1:nrow(matched_pooled_noag)
    matched_pooled_nlo$rowid <- 1:nrow(matched_pooled_nlo)
    
###############################################
# Descriptive Figures and Tables ##############
###############################################
    
# scatter plot
  
  cairo_ps("./results/region_18th_scatter.eps",width=6,height=3.5,fallback_resolution = 1200)
  pooled_scatter <- ggplot(pooled_ag,aes(x=score_mean,y=vote_prohibition,group=Region,colour=Region))+
    geom_jitter(height=0.03,width=0,size=2,alpha=0.2)+
    geom_smooth(se=FALSE,size=3,method="loess")+
    geom_vline(aes(xintercept=0.5),size=2,linetype=2,colour="black")+
    ylab("Vote on Prohibition")+
    xlab("District Preference")+
    scale_color_manual(values=c("black","grey50"))+
    scale_fill_manual(values=c("black","grey50"))+
    theme_minimal()
  print(pooled_scatter)
  dev.off()
  
# histogram of matched preferences
  
  cairo_ps("./results/matched_preference_histogram.eps",width=3.25,height=2,fallback_resolution = 1200)
  pref_hist <- 
    ggplot(data=matched_pooled,aes(x=score_mean,group=Region,fill=Region,colour=Region))+
    geom_density(size=0.75,alpha=0.3)+
    xlim(0,1)+
    xlab("Share District Pro-Prohibition")+
    ylab("Density")+
    scale_color_manual(values=c("black","grey50"))+
    scale_fill_manual(values=c("black","grey50"))+
    theme_minimal()
  print(pref_hist)
  dev.off()
  
# histogram of non-matched preferences
  
  cairo_ps("./results/sl_preference_histogram.eps",width=3.25,height=2,fallback_resolution = 1200)
  pref_hist <- 
    ggplot(data=pooled_ag,aes(x=score_mean,group=Region,fill=Region,colour=Region))+
    geom_density(size=0.75,alpha=0.3)+
    xlim(0,1)+
    xlab("Share District Pro-Prohibition")+
    ylab("Density")+
    scale_color_manual(values=c("black","grey50"))+
    scale_fill_manual(values=c("black","grey50"))+
    theme_minimal()
  print(pref_hist)
  dev.off()

# density plot of vote on 18th amendment
  
  cairo_ps("./results/sl_vote_histogram.eps",width=6,height=3.5,fallback_resolution = 1200)
  pref_hist <- 
    ggplot(data=pooled_ag,aes(x=vote_prohibition,group=Region,fill=Region,colour=Region))+
    geom_density(size=0.75,alpha=0.3)+
    xlim(0,1)+
    xlab("Share Votes Pro-Prohibition")+
    ylab("Density")+
    scale_color_manual(values=c("black","grey50"))+
    scale_fill_manual(values=c("black","grey50"))+
    theme_minimal()
  print(pref_hist)
  dev.off()
  
# summary statistics table
  
  sumstat_df <- as.data.frame(pooled_ag[,c("vote_prohibition","score_mean","n_votes","south")])
  
  sumstats_sg <- stargazer(sumstat_df,summary=TRUE,
                           title="Summary Statistics, 18th Amendment Ratification",
                           label="sumstats_amendment",
                           covariate.labels = c("Prohibition Roll Call Score",
                                                "District Preference on Prohibition",
                                                "Number of Prohibition Votes",
                                                "South"),
                           summary.stat=c("mean","median","sd","min","max","n"),
                           out="./results/sumstats_amendment.tex")
  
# matched summary statistics table
  
  sumstat_matched_df <- as.data.frame(matched_pooled[,c("vote_prohibition","score_mean","n_votes","south")])
  
  sumstats_matched_sg <- stargazer(sumstat_matched_df,summary=TRUE,
                           title="Summary Statistics, 18th Amendment Ratification, Matched Sample",
                           label="sumstats_amendment_matched",
                           covariate.labels = c("Prohibition Roll Call Score",
                                                "District Preference on Prohibition",
                                                "Number of Prohibition Votes",
                                                "South"),
                           summary.stat=c("mean","median","sd","min","max","n"),
                           out="./results/sumstats_amendment_matched.tex")
  
##############################################
# Analysis ###################################
##############################################
  
# regressions
  
  nocov18 <- felm(vote_prohibition~score_mean*south|0|0|rowid,data=pooled_ag)
  
  nocov18_weight <- felm(vote_prohibition~score_mean*south|0|0|rowid,weight=pooled_ag$n_votes,data=pooled_ag)
  
  nocov18_matched <- felm(temp_outcome~score_mean*south|0|0|rowid,data=matched_pooled)
  
  nocov18_weight_matched <- felm(temp_outcome~score_mean*south|0|0|rowid,weight=matched_pooled$n_votes,data=matched_pooled)
  
  amendment_sg <- stargazer(nocov18,nocov18_weight,nocov18_matched,nocov18_weight_matched,
                            notes.append = FALSE,notes.label = "",
                            notes="\\parbox[t]{0.65\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                        robust standard errors in parentheses. Weighted models are weighted according to the number of votes 
                        cast by each district's representatives. $^{*}$p$<$0.05 (two-tailed test).}",
                            add.lines = list(c("Weights","","\\checkmark","","\\checkmark")),
                            keep.stat = c("adj.rsq","n"),
                            star.char=c("*"),star.cutoffs = c(0.05),digits=3,digits.extra=0,
                            table.placement = "!ht",
                            dep.var.caption = "Share Votes Pro-Prohibition",
                            dep.var.labels = c("Unmatched","Matched"),
                            covariate.labels = c("District Preference",
                                                 "South","Preference $\\times$ South",
                                                 "Constant"),
                            label="amendment",no.space=T,
                            table.layout ="-ld-#-t-as-n",
                            title="Constituent Preferences and 18th Amendment Ratification")
  
  cat(amendment_sg, sep = '\n', file = "./results/amendment.tex")
  
# regressions, state fixed effects
  
  nocov18_fe <- felm(vote_prohibition~score_mean*south-south|state|0|rowid,data=pooled_ag)
  
  nocov18_fe_weight <- felm(vote_prohibition~score_mean*south-south|state|0|rowid,weight=pooled_ag$n_votes,data=pooled_ag)
  
  nocov18_fe_matched <- felm(temp_outcome~score_mean*south-south|state|0|rowid,data=matched_pooled)
  
  nocov18_fe_weight_matched <- felm(temp_outcome~score_mean*south-south|state|0|rowid,weight=matched_pooled$n_votes,data=matched_pooled)
  
  amendment_fe_sg <- stargazer(nocov18_fe,nocov18_fe_weight,nocov18_fe_matched,nocov18_fe_weight_matched,
                            notes.append = FALSE,notes.label = "",
                            notes="\\parbox[t]{0.6\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                            robust standard errors in parentheses. Weighted models are weighted according to the number of votes 
                            cast by each district's representatives. $^{*}$p$<$0.05 (two-tailed test).}",
                            add.lines = list(c("Weights","","\\checkmark","","\\checkmark"),
                                             c("State FE","\\checkmark","\\checkmark","\\checkmark","\\checkmark")),
                            keep.stat = c("adj.rsq","n"),
                            star.char=c("*"),star.cutoffs = c(0.05),digits=3,digits.extra=0,
                            table.placement = "!ht",
                            dep.var.caption = "Share Votes Pro-Prohibition",
                            dep.var.labels = c("Unmatched","Matched"),
                           covariate.labels = c("District Preference",
                                                "Preference $\\times$ South"),
                            label="amendment_fe",
                            table.layout ="-ld-#-t-as-n",
                            title="18th Amendment Ratification, State Fixed Effects")
  
  cat(amendment_fe_sg, sep = '\n', file = "./results/amendment_fe.tex")
  
# regressions, no local option
  
  nocov18_nlo <- felm(vote_prohibition~score_mean_nlo*south|0|0|rowid,data=pooled_ag)
  
  nocov18_weight_nlo <- felm(vote_prohibition~score_mean_nlo*south|0|0|rowid,weight=pooled_ag$n_votes,data=pooled_ag)
  
  nocov18_matched_nlo <- felm(temp_outcome~score_mean_nlo*south|0|0|rowid,data=matched_pooled_nlo)

  nocov18_weight_matched_nlo <- felm(temp_outcome~score_mean_nlo*south|0|0|rowid,weight=matched_pooled_nlo$n_votes,data=matched_pooled_nlo)
  
  amendment_nlo_sg <- stargazer(nocov18_nlo,nocov18_weight_nlo,nocov18_matched_nlo,nocov18_weight_matched_nlo,
                            notes.append = FALSE,notes.label = "",
                            notes="\\parbox[t]{0.65\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                        robust standard errors in parentheses. Weighted models are weighted according to the number of votes 
                        cast by each district's representatives. $^{*}$p$<$0.05 (two-tailed test).}",
                            add.lines = list(c("Weights","","\\checkmark","","\\checkmark")),
                            keep.stat = c("adj.rsq","n"),
                            star.char=c("*"),star.cutoffs = c(0.05),digits=3,digits.extra=0,
                            table.placement = "!ht",
                            dep.var.caption = "Share Votes Pro-Prohibition",
                            dep.var.labels = c("Unmatched","Matched"),
                            covariate.labels = c("District Preference", 
                                                 "South","Preference $\\times$ South",
                                                 "Constant"),
                            label="amendment_nlo",
                            table.layout ="-ld-#-t-as-n",
                            title="18th Amendment Ratification, No Local Option")
  
  cat(amendment_nlo_sg, sep = '\n', file = "./results/amendment_nlo.tex")
  
# drop counties with too many votes
  
  # share of counties with only one vote
  
    nrow(pooled_ag[pooled_ag$n_votes==1,])/nrow(pooled_ag)
  
  onevote <- felm(vote_prohibition~score_mean*south|0|0|rowid,data=pooled_ag[pooled_ag$n_votes==1,])
  
  onevote_matched <- felm(temp_outcome~score_mean*south|0|0|rowid,data=matched_pooled[matched_pooled$n_votes==1,])
  
  onevote_sg <- stargazer(onevote,onevote_matched,
                            notes.append = FALSE,notes.label = "",
                            notes="\\parbox[t]{0.55\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                            robust standard errors in parentheses. Matching based on full sample. $^{*}$p$<$0.05 (two-tailed test).}",
                            add.lines = list(c("Weights","","\\checkmark","","\\checkmark")),
                            keep.stat = c("adj.rsq","n"),
                            star.char=c("*"),star.cutoffs = c(0.05),digits=3,digits.extra=0,
                            table.placement = "!ht",
                            dep.var.caption = "Share Votes Pro-Prohibition",#font.size="footnotesize",
                            no.space=T,
                            dep.var.labels = c("Unmatched","Matched"),
                            covariate.labels = c("District Preference", 
                                                 "South","Preference $\\times$ South",
                                                 "Constant"),
                            label="onevote",
                            table.layout ="-ld-#-t-as-n",
                            title="18th Amendment Ratification, One-Vote County Sample")
  
  cat(onevote_sg, sep = '\n', file = "./results/onevote.tex")

# leave one state out  
  
  pooled_ag <- dplyr::arrange(pooled_ag,desc(south),state,county)
  
  state_out <- list()
  
  for(j in 1:length(unique(pooled_ag$state))){
    
    bit <- pooled_ag[pooled_ag$state!=sort(unique(pooled_ag$state))[j],]
    
    reg   <- felm(vote_prohibition~score_mean*south|0|0|rowid,data=bit)
    
    state_out[[j]] <- c(sort(unique(as.character(pooled_ag$state)))[j],coef(reg)["score_mean:south"],
                        sqrt(reg$clustervcv["score_mean:south","score_mean:south"]))
    
  }
  
  state_out <- as.data.frame(do.call(rbind,state_out))
  colnames(state_out) <- c("state","coef","se")
  state_out$coef <- as.numeric(as.character(state_out$coef))
  state_out$se <- as.numeric(as.character(state_out$se))
  state_out$south <- 0; state_out[is.element(state_out$state,
                                             c("AR","AL","GA","FL","NC","SC","TX","VA","TN","LA","MS")),"south"] <- 1
  state_out <- dplyr::arrange(state_out,desc(south),state)
  state_out$Region <- factor(ifelse(state_out$south==1,"South","North"))
  state_out$state <- factor(state_out$state,
                            levels=c("AZ","AR","CA","CO",
                                    "FL","ID","IN","IA",
                                    "MI","MT","NE","ND",
                                    "OH","OK","OR","SC",
                                    "SD","TX","UT","VA","WA","WV"))
  
  
  cairo_ps("./results/omitted_state_amendment.eps",width=7,height=3.5,fallback_resolution = 1200)
  omit_state_plot <- ggplot(data=state_out,aes(x=state,y=coef,colour=Region))+
    geom_point(size=2)+
    geom_errorbar(aes(ymax=coef+1.96*se,ymin=coef-1.96*se),size=1.25)+
    scale_color_manual(values=c("black","grey70"))+
    scale_fill_manual(values=c("black","grey70"))+
    geom_hline(yintercept=coef(nocov18)["score_mean:south"],
              size=2,colour="indianred",linetype=2)+
    ylab("Coefficient on 'South X District Preference'")+
    xlab("Omitted State")+
    theme_minimal()
  print(omit_state_plot)
  dev.off()
  
# logit
  
  pooled$temp_outcome <- pooled$vote_prohibition
  
  ols <- glm(vote_prohibition~score_mean*south,family=gaussian,data=pooled)
  ols_vcov <- cluster.vcov(ols,pooled$county)
  
  ols_match <- glm(vote_prohibition~score_mean*south,family=gaussian,data=matched_pooled_noag)
  ols_match_vcov <- cluster.vcov(ols_match,matched_pooled_noag$county)
  
  ols_maj <- glm(vote_prohibition~I(score_mean>=0.5)*south,family=gaussian,data=pooled)
  ols_maj_vcov <- cluster.vcov(ols_maj,pooled$county)
  
  ols_maj_match <- glm(vote_prohibition~I(score_mean>=0.5)*south,family=gaussian,data=matched_pooled_noag)
  ols_maj_match_vcov <- cluster.vcov(ols_maj_match,matched_pooled_noag$county)
  
  logit <- glm(temp_outcome~score_mean*south,family=binomial,data=pooled)
  logit_vcov <- cluster.vcov(logit,pooled$county)
  
  logit_match <- glm(temp_outcome~score_mean*south,family=binomial,data=matched_pooled_noag)
  logit_match_vcov <- cluster.vcov(logit_match,matched_pooled_noag$county)
  
  logit_maj <- glm(temp_outcome~I(score_mean>=0.5)*south,family=binomial,data=pooled)
  logit_maj_vcov <- cluster.vcov(logit_maj,pooled$county)
  
  logit_maj_match <- glm(temp_outcome~I(score_mean>=0.5)*south,family=binomial,data=matched_pooled_noag)
  logit_maj_match_vcov <- cluster.vcov(logit_maj_match,matched_pooled_noag$county)
  
  logit_sg <- stargazer(ols,ols_match,ols_maj,ols_maj_match,logit,logit_match,logit_maj,logit_maj_match,
                         se=list(sqrt(diag(ols_vcov)),sqrt(diag(ols_match_vcov)),sqrt(diag(ols_maj_vcov)),sqrt(diag(ols_maj_match_vcov)),
                                 sqrt(diag(logit_vcov)),sqrt(diag(logit_match_vcov)),sqrt(diag(logit_maj_vcov)),sqrt(diag(logit_maj_match_vcov))),
                     notes.append = FALSE,notes.label = "",
                     notes="\\parbox[t]{0.995\\textwidth}{\\footnotesize \\textit{Note}: Entries are coefficients with 
                             district-clustered standard errors in parentheses. Models 1-4 are linear probability models, and Models 5-8 are logit models.
                             $^{*}$p$<$0.05 (two-tailed test).}",
                     add.lines = list(c("Matched Sample","","\\checkmark","","\\checkmark","","\\checkmark","","\\checkmark")),
                     keep.stat = c("n"),font.size="footnotesize",column.sep.width = "0pt",
                     star.char=c("*"),star.cutoffs = c(0.05),digits=3,digits.extra=0,no.space=T,
                     table.placement = "!ht",
                     dep.var.caption = "Vote Pro-Prohibition",
                     dep.var.labels = c("Linear Probability Model","Logit"),
                      covariate.labels = c("District Preference",
                                           "District Maj. for Prohibition",
                                           "South",
                                           "Preference $\\times$ South",
                                           "Maj. for Prohibition $\\times$ South",
                                           "Constant"),
                     label="logit",
                     table.layout ="-ld-#-t-as-n",
                     title="18th Amendment Ratification, Legislator-Level Models")
  
  cat(logit_sg, sep = '\n', file = "./results/logit.tex")

  
# majority regressions
  
  nocov18_maj <- felm(vote_prohibition~I(score_mean>=0.5)*south|0|0|rowid,data=pooled_ag)
  
  nocov18_weight_maj <- felm(vote_prohibition~I(score_mean>=0.5)*south|0|0|rowid,
                         weight=pooled_ag$n_votes,data=pooled_ag)
  
  nocov18_matched_maj <- felm(temp_outcome~I(score_mean>=0.5)*south|0|0|rowid,data=matched_pooled)
  
  nocov18_weight_matched_maj <- felm(temp_outcome~I(score_mean>=0.5)*south|0|0|rowid,
                                 weight=matched_pooled$n_votes,data=matched_pooled)
  
  amendment_maj_sg <- stargazer(nocov18_maj,nocov18_weight_maj,nocov18_matched_maj,nocov18_weight_matched_maj,
                            notes.append = FALSE,notes.label = "",
                            notes="\\parbox[t]{0.85\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                            robust standard errors in parentheses. Sample is all districts with district preference measures less than 0.45
                            or greater than 0.55. Weighted models are weighted according to the number of votes 
                            cast by each district's representatives. $^{*}$p$<$0.05 (two-tailed test).}",
                            add.lines = list(c("Weights","","\\checkmark","","\\checkmark")),
                            keep.stat = c("adj.rsq","n"),
                            star.char=c("*"),star.cutoffs = c(0.05),digits=3,digits.extra=0,
                            table.placement = "!ht",
                            dep.var.caption = "Share Votes Pro-Prohibition",
                            no.space=TRUE,
                            dep.var.labels = c("Unmatched","Matched"),
                            covariate.labels = c("District Maj. for Prohibition",
                                                "South","District Maj. for Prohibition $\\times$ South",
                                                "Constant"),
                            label="amendment_maj",
                            table.layout ="-ld-#-t-as-n",
                            title="Majority Preferences and 18th Amendment Ratification")
  
  cat(amendment_maj_sg, sep = '\n', file = "./results/amendment_maj.tex")

# lopping off the middle
  
  lop_sequence <- seq(from=0,to=0.3,by=0.025)
  
  lop_coef <- list()
  
  for(j in 1:length(lop_sequence)){
    
    loop_weight_reg <- felm(vote_prohibition~score_mean*south|0|0|rowid,
                            data=pooled_ag[abs(pooled_ag$score_mean-0.5)>lop_sequence[j],])
    
    lop_coef[[j]] <- summary(loop_weight_reg)$coefficients["score_mean:south",]
    
  }
  
  lop_coef <- as.data.frame(cbind(lop_sequence,do.call(rbind,lop_coef)))
  
  
  cairo_ps("./results/sl_pref_lopped_estimates.eps",width=6,height=3.5,fallback_resolution = 1200)
  lopped_plot <- ggplot(lop_coef,aes(x=lop_sequence,y=Estimate))+
    geom_point(size=2)+
    geom_errorbar(aes(ymin=Estimate-1.96*`Cluster s.e.`,ymax=Estimate+1.96*`Cluster s.e.`),width=0.015)+
    geom_hline(yintercept=0,linetype=2,colour="indianred",size=2)+
    xlab("Lopped Distance")+
    ylab("Coefficient on 'District Preference X South'")+
    theme_minimal()
  print(lopped_plot)
  dev.off()
  
  # create a "solid north" variable
  
  dubin <- read.dta("./seat_numbers.dta")
  dubin_early <- dubin[dubin$year<=1920 & dubin$year>=1910,]
  dubin_early <- as.data.frame(dubin_early %>% group_by(state) %>% summarise(n_dem = sum(lhd+uhd,na.rm=TRUE),
                                                                             n_rep = sum(lhr+uhr,na.rm=TRUE),
                                                                             n_tot = sum(lht+uht,na.rm=TRUE)))
  dubin_early$rep_share <- dubin_early$n_rep/dubin_early$n_tot
  dubin_early$dem_share <- dubin_early$n_dem/dubin_early$n_tot
  dubin_early$south <- ifelse(is.element(dubin_early$state,c("AL","AR","FL","GA","LA","MS","NC","SC","TN","TX","VA")),1,0)
  
  pooled_ag$solid_north <- ifelse(is.element(pooled_ag$state,dubin_early$state[(dubin_early$rep_share>=0.7 | 
                                                                                 dubin_early$dem_share>=0.7) & dubin_early$south==0]),1,0)
  
  nocov18_north <- felm(vote_prohibition~score_mean*(south+solid_north)|0|0|rowid,
                      data=pooled_ag)
  
  nocov18_weight_north <- felm(vote_prohibition~score_mean*(south+solid_north)|0|0|rowid,
                             weight=pooled_ag$n_votes,data=pooled_ag)
  
  amendment_north_sg <- stargazer(nocov18_north,nocov18_weight_north,
                                notes.append = FALSE,notes.label = "",
                                notes="\\parbox[t]{0.675\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                                robust standard errors in parentheses. `Solid North` is an indicator for whether 70\\% or more of a state's state legislators
                                were of one party between 1910 and 1920.
                                Weighted models are weighted according to the number of votes 
                                cast by each district's representatives. $^{*}$p$<$0.05 (two-tailed test).}",
                                add.lines = list(c("Weights","","\\checkmark","","\\checkmark")),
                                keep.stat = c("adj.rsq","n"),
                                star.char=c("*"),star.cutoffs = c(0.05),digits=3,digits.extra=0,
                                table.placement = "!ht",no.space=T,
                                dep.var.labels = "Share Votes Pro-Prohibition",
                                covariate.labels = c("District Preference",
                                                   "South","'Solid North'","District Preference $\\times$ South",
                                                   "District Preference $\\times$ 'Solid North'",
                                                   "Constant"),
                                label="sn_sl",
                                table.layout ="-ld-#-t-as-n",
                                title="18th Amendment Ratification, `Solid North'")
  
  cat(amendment_north_sg, sep = '\n', file = "./results/amendment_north.tex")
  
# terciles  
  
  amendment_terciles   <- felm(vote_prohibition~terciles*south-1-south|0|0|rowid,data=pooled_ag)
  
  amendment_weights_terciles <- felm(vote_prohibition~terciles*south-1-south|0|0|rowid,data=pooled_ag,weights=pooled_ag$n_votes)
  
  amendment_matched_terciles   <- felm(temp_outcome~terciles*south-1-south|0|0|rowid,data=matched_pooled)
  
  amendment_weights_matched_terciles <- felm(temp_outcome~terciles*south-1-south|0|0|rowid,data=matched_pooled,weights=matched_pooled$n_votes)
  
  amendment_tercile_sg <- stargazer(amendment_terciles,amendment_weights_terciles,amendment_matched_terciles,amendment_weights_matched_terciles,
                                    notes.append = FALSE,notes.label = "",
                                    notes=paste("\\parbox[t]{0.7\\textwidth}{\\footnotesize \\textit{Note}: Entries are linear regression coefficients with 
                                                robust standard errors in parentheses. Models estimated without omitted category. Weighted models are weighted according to the number of votes 
                                                cast by each district's representatives. Continuous measure of district preference has been binned into 
                                                terciles based on the full sample, pooled across South and non-South. First tercile is between ", round(quantile(pooled_ag$score_mean,0)[1],2),
                                                " and ", round(quantile(pooled_ag$score_mean,0.33333333)[1],2),", the second is between ", round(quantile(pooled_ag$score_mean,0.33333333)[1],2), " and ", 
                                                round(quantile(pooled_ag$score_mean,0.66666666)[1],2),", and the third between ", round(quantile(pooled_ag$score_mean,0.33333333)[1],2), " and ",
                                                round(quantile(pooled_ag$score_mean,1)[1],2),". ``Joint Test $p$'' is the $p$ value for a joint $F$ test for the hypothesis that ``South $\\times$ Tercile One,'' ``South $\\times$ Tercile Two,''
                                                and ``South $\\times$ Tercile Three'' are each zero. $^{*}p<$0.05 (two-tailed test).}",sep=""),
                                    add.lines = list(c("Weights","","\\checkmark","","\\checkmark"),
                                                     c("Joint Test $p$",
                                                       format(round(waldtest(amendment_terciles,~`terciles1:south`|`terciles2:south`|`terciles3:south`)["p.F"],2),nsmall=2),
                                                       format(round(waldtest(amendment_weights_terciles,~`terciles1:south`|`terciles2:south`|`terciles3:south`)["p.F"],2),nsmall=2),
                                                       format(round(waldtest(amendment_matched_terciles,~`terciles1:south`|`terciles2:south`|`terciles3:south`)["p.F"],2),nsmall=2),
                                                       format(round(waldtest(amendment_weights_matched_terciles,~`terciles1:south`|`terciles2:south`|`terciles3:south`)["p.F"],2),nsmall=2))),
                                    keep.stat = c("adj.rsq","n"),
                                    star.char=c("*"),star.cutoffs = c(0.05),digits=3,digits.extra=0,no.space = T,
                                    table.placement = "!ht",
                                    dep.var.caption = "Share Votes Pro-Prohibition",
                                    dep.var.labels = c("Unmatched","Matched"),
                                    covariate.labels = c("Tercile One","Tercile Two","Tercile Three",
                                                         "South $\\times$ Tercile One","South $\\times$ Tercile Two","South $\\times$ Tercile Three"),
                                    label="amendment_tercile",
                                    table.layout ="-ld-#-t-as-n",
                                    title="18th Amendment Ratification, Binned Preferences")
  
  cat(amendment_tercile_sg, sep = '\n', file = "./results/amendment_tercile.tex")
  
  
  
  
